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Signal-Selective  Time-Dilference-of-Arrival 
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Sources  in  Highly  Cormptive  Environments,  Part  I: 

Theory  and  Method 
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Abstract— For  the  problem  of  estimating  time  difference  of 
arrival  (TDOA)  of  radio  waves  impinging  on  a  pair  of  antennas 
for  the  purpose  of  passively  locating  the  source  of  a  communi¬ 
cations  or  telemetry  signal  in  the  presence  of  interfering  signals 
and  noise,  a  new  class  of  signal-selective  algorithms  that  is 
highly  tolerant  to  interference  and  noise  is  introduced.  In  part 
I  of  this  two-part  paper,  the  background  theory  of  cyclosta¬ 
tionary  signals  is  presented  and  applied  to  the  design  of  various 
new  TDOA  methods.  In  part  II,  algorithmic  implementations 
are  described  and  their  performance  capabilities  are  assessed 
by  analysis  and  simulation.  By  virtue  of  the  fact  that  the  mul¬ 
tiple-signal  resolution  problem  is  essentially  eliminated  in  many 
applications  by  the  signal  selectivity  of  the  algorithms,  two  per¬ 
formance  advantages  are  gained:  I)  the  practicality  of  source 
location  based  on  time-difference  measurements  with  relatively 
closely  spaced  antennas  is  substantially  enhanced,  and  2)  the 
problem  of  sorting  through  multiple  TDOA  estimates,  resulting 
from  multiple  interferers,  for  the  estimate  corresponding  to  a 
particular  signal  of  interest  is  eliminated.  These  new  algo¬ 
rithms  exhibit  their  signal  selectivity  regardless  of  the  extent  of 
temporal,  spectral,  or  spatial  overlap  among  received  signals. 
It  is  only  required  that  the  signal  of  interest  have  a  known  (or 
measurable)  analog  carrier  frequency  or  digital  keying  rate  that 
is  distinct  from  those  of  all  interfering  signals.  Yet  the  com¬ 
putational  complexity  of  these  algorithms  is  comparable  to  that 
of  conventional  generalized  cross-correlation  algorithms. 


I.  Introduction 

A.  Background 

The  problem  of  passively  locating  the  source  of  a  prop¬ 
agating  electromagnetic  wave  consisting  of  a  man-made 
signal  has  a  number  of  applications  including  direction 
finding  for  navigation  by  land,  air,  or  sea,  tracking  of 
moving  emitters  for  surveillance  (e.g.,  for  research,  law 
enforcement,  or  national  security),  monitoring  and  locat- 
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ing  illegal  and/or  hostile  communicators  (e.g.,  contra¬ 
band,  violators  of  communications  regulations,  or  ene¬ 
mies  in  a  battlefield),  military  reconnaissance  and 
intelligence,  and  so  on.  Unlike  active  source-location  sys¬ 
tems,  like  radar,  sonar,  and  lidar,  which  transmit  signals 
and  then  process  the  received  reflections  from  the  objects 
to  be  located,  passive  systems  simply  process  whatever 
signals  are  emitted  from  the  object  to  be  located. 

In  many  of  these  passive  source-location  applications, 
the  receivers  in  the  source-location  system  are  subject  to 
a  variety  of  types  of  electromagnetic  interference  and 
noise,  including  natural  and  manmade  signals  other  than 
the  signal  of  interest  (SOI).  These  signals  not  of  interest 
(SNOI)  can  severely  degrade  the  performance  of  conven¬ 
tional  source-location  systems  when  they  are  present  at 
the  same  time  and  also  occupy  the  same  spectral  band  as 
the  SOI.  This  can  be  particularly  problematic  when  the 
SOI  is  weak  relative  to  the  SNOI  and/or  noise.  The  prob¬ 
lem  is  further  exacerbated  when  the  locations  of  the 
sources  of  the  SNOI  are  unknown  and/or  close  to  that  of 
the  SOI. 

The  purpose  of  this  paper  is  to  present  a  new  class  of 
passive  source-location  methods  that  exploit  an  inherent 
signal  property,  called  cyclostationarity,  that  is  charac¬ 
teristic  of  man-made  signals  used  for  communications  and 
telemetry  to  obtain  substantial  tolerance  to  all  types  of 
interference  and  noise  (except  for  some  interfering  signals 
that  are  intentionally  designed  to  exhibit  the  same  cyclo- 
stationarity,  e.g.,  in  communication  networks).  Like  most 
conventional  methods  (cf.  [2])  that  require  some  degree 
of  tolerance  to  noise,  the  new  methods  are  based  on  cross 
correlation  of  time-shifted  measurements  of  data  from 
multiple  receivers.  However,  unlike  conventional  meth¬ 
ods,'  the  new  methods  cross  correlate  frequency-shifted 

'Although  cross-ambiguity  methods  (cf.  [3]),  which  compensate  for  fre¬ 
quency  difference  of  arrival  due  to  Doppler  effects,  do  cross  correlate  fre¬ 
quency-shifted  data,  the  frequency  shifts  are  relatively  small  since  they 
correspond  to  Doppler  shifts,  whereas  the  frequency  shafts  used  in  the  new 
methods  are  much  larger  than  typical  Doppler  shifts  a-.d  accomplish  an 
entirely  different  task. 
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as  well  as  time-shifted  versions  of  the  received  data  in 
order  to  exploit  the  unique  cyclostationarity  property  of 
the  SOL. 

Within  the  general  class  of  source-location  methods  that 
are  based  on  cross-correlation  measurements,  there  are 
two  distinct  subclasses:  There  are  those  methods  that  are 
designed  for  an  array  of  closely  spaced  antenna  elements 
on  a  single  platform  (e.g.,  a  ship,  aircraft,  satellite, 
ground  vehicle,  or  fixed  ground  station),  for  which  the 
element  spacing  is  typically  less  than  half  a  wavelength 
for  all  signals  received  and  phase-alignment  methods  for 
beam/null  steering  are  employed;  and  there  are  those 
methods  that  are  designed  for  two  or  three  widely  spaced 
antenna  elements,  each  often  (but  not  always)  on  a  sepa¬ 
rate  platform,  where  time-difference  measurements  are 
used  to  obtain  location  information.  Whereas  the  former 
class  of  methods  (often  called  direction  finding  methods) 
exploits  phase  differences  (less  than  t  rad  in  order  to  avoid 
ambiguities),  from  element  to  element,  of  relatively  nar¬ 
row-band  signals  (or  wide-band  signals  decomposed  into 
narrow-band  components)  to  estimate  angle  of  arrival,  the 
latter  class  of  methods  (often  called  time-difference  lo¬ 
cation  methods)  exploits  relative  time  differences  (the 
larger  the  better  since  root-mean-squared  error  of  location 
is  approximately  inversely  proportional  to  the  separation 
between  platforms)  of  preferably  wide-band  signals  from 
platform  to  platform  to  estimate  location  (both  angle  of 
arrival  and  range).  In  practice,  however,  both  methods 
can  be  applied  to  all  bandwidths  used  in  typical  commu¬ 
nications  and  telemetry  systems. 

The  requirements  on  accuracy  and  spatial  resolution  ca¬ 
pabilities  of  array-based  methods  become  more  stringent 
as  the  distance  between  each  source  to  be  located  and  the 
reception  platform  increases,  since  this  decreases  differ¬ 
ences  between  angles  of  arrival.  In  contrast,  the  require¬ 
ments  on  accuracy  and  temporal  resolution  capabilities  of 
time-difference-of-arrival  (TDOA)-based  methods  be¬ 
come  less  stringent  as  the  separation  between  reception 
platforms  increases,  since  this  increases  differences  be¬ 
tween  times  of  arrival. 

The  need  for  high  resolution  arises  primarily  when  (rel¬ 
atively)  closely  spaced  multiple  sources  give  rise  to  mul¬ 
tiple  received  signals  that  cannot  be  separated  by  prepro¬ 
cessing  (prior  to  processing  for  location).  For  instance, 
TDOA’s  of  multiple  signals  that  are  not  separated  by  more 
than  the  widths  of  their  cross-correlation  peaks  (whose 
locations  on  the  time-delay  axis  correspond  to  the 
TDOA’s)  usually  cannot  be  resolved  by  conventional 
TDOA-based  methods.  To  minimize  this  problem,  the 
distance  between  platforms  is  typically  made  as  large  as 
is  practically  feasible  so  that  the  magnitudes  of  the 
TDOA’s  (which  are  proportional  to  the  distance  between 
platforms)  will  be  as  large  as  possible,  thereby  minimiz¬ 
ing  the  overlap  of  adjacent  peaks  (whose  widths  dejTend 
only  on  the  signal  bandwidth— being  inversely  propor¬ 
tional— not  on  the  distance  between  platforms).  Specifi¬ 
cally,  in  planar  space  the  approximate  minimum  resolv¬ 


able  separation  A6  in  angles  of  arrival  (AOA)  of  two 
signals  in  the  far  field  of  the  platform  pair  with  approxi¬ 
mately  linear  wavefronts  can  be  derived  as  follows; 

^  <  (D|  -  Dt  I  =  -  (sin  d\  -  sin  6-,  | 

B  c 


where  D,  and  D2  are  the  two  TDOA’s  in  seconds,  and 
62  are  the  corresponding  AOA’s  in  radians  relative  to 
broadside,  L  is  the  distance  (in  meters)  between  receivers, 
B  is  the  signal  bandwidth  (in  Hertz),  and  c  is  the  speed  of 
propagation  (3  x  10*  m/s).  For  Afl  4  6,  -  62  small,  we 
can  use  sin  (A6/2)  =  Ad/2  in  this  inequality  to  obtain 


where 


-  LBlcos  d\ 


0  4  (0,  +  62)/2. 


For  example,  for  L  =  3000  m,  B  =  1  MHz,  and  6  =  60 
degrees,  we  obtain  A6^i„  =  12  degrees.  This  relatively 
poor  resolution  becomes  worse  as  either  L  or  B  is  de¬ 
creased. 

The  best  performing  array-based  methods  attempt  to 
circumvent  this  resolution  problem  in  locating  multiple 
signal  sources  by  simultaneously  estimating  multiple  an¬ 
gles  of  arrival  rather  than  individually  estimating  the  an¬ 
gle  of  arrival  of  each  signal  (as  commonly  done  by  con¬ 
ventional  beamformers  and  TDOA-based  methods)  [7]. 
Although  the  array-based  methods  offer  the  advantage  of 
high  spatial  resolution  (with  respect  to  the  spatial  extent 
of  the  source-location  system— the  array  size),  and  the 
ability  to  simultaneously  locate  a  number  of  signals  up  to 
one  less  than  the  number  of  elements  in  the  array,  their 
complexity  is  typically  much  higher  than  that  of  the  time- 
difference-of-arrival-based  methods  because  of  the  need 
for  measurement,  storge,  and  usage  of  large  amounts  of 
array  calibration  data  (the  recently  proposed  ESPRIT 
method  is  an  exception  because  of  its  special  array  design 
[4]),  and  because  of  the  use  of  computationally  intensive 
algorithms:  the  best  performing  algorithms  require  sin¬ 
gular  value  decomposition  of  cross-correlation  matrices 
of  possibly  high  dimension,  and/or  require  the  solution  of 
a  multidimensional  optimization  problem  [8],  [9]. 


B.  Overview 

The  new  cyclostationarity-exploiting  TDOA-based 
methods  presented  in  this  paper  significantly  alter  this 
tradeoff  between  highly  sophisticated  high-resolution  ar¬ 
ray-based  methods  and  the  relatively  simple  TDOA-based 
methods  that  require  widely  separated  multiple  platforms 
by  eliminating  the  resolution-limitation  problem  for 
TDOA.  Thpt  is.  because  the  spectral  correlations  used 
make  the  new  methods  signal  selective  (in  spite  of  tern- 
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poral  and  spectral  overlap  between  the  SOI  and  SNOI), 
they  usually  do  not  need  to  resolve  multiple  crosscorre¬ 
lation  peaks  along  the  time-delay  axis.  As  a  result  of  this 
novel  signal  selectivity,  the  requirements  on  separation  of 
platforms  can  be  eased  considerably.  In  fact,  single-plat¬ 
form  TDOA-based  source-location  systems  in  environ¬ 
ments  with  SNOI  become  more  practically  feasible  with 
the  new  methods.  Also,  the  problem  of  deciding  which 
peak  among  multiple  cross-correlation  peaks  is  associated 
with  the  SOI,  and  which  with  the  SNOI,  is  eliminated. 
The  only  requirement  of  the  new  methods  is  that  they  must 
know  a  carrier  frequency,  or  keying  rate,  or  possibly  some 
other  cycle  frequency  associated  with  the  cyclostationar- 
ity  of  the  random  SOI,  although  such  cycle  frequencies 
can  be  estimated  using  the  same  data  to  be  used  for  TDOA 
estimation.  These  cycle  frequencies  are  the  frequency 
shifts  used  to  obtain  signal  selectivity. 

The  objectives  of  this  two-part  paper  are  threefold;  The 
first  objective,  which  is  met  in  part  I,  is  to  introduce  a 
variety  of  new  methods  for  signal-selective  TDOA  esti¬ 
mation;  the  second  objective,  which  is  met  in  part  II  [1], 
is  to  present  the  first  quantitative  results  on  a  comparative 
performance  study  of  some  of  these  algorithms  based  on 
simulations.  The  results  show  that  although  all  the  algo¬ 
rithms  tested  do  indeed  offer  the  advantages  described  in 
this  introduction,  which  accrue  with  signal  selectivity, 
some  algorithms  are  superior  in  the  sense  of  yielding  equal 
mean  squared  error  (MSE)  in  the  TDOA  estimates  with 
substantially  less  received  data.  In  fact,  some  of  these  al¬ 
gorithms  similarly  outperform  conventional  TDOA  algo¬ 
rithms  even  in  the  absence  of  SNOI,  that  is,  when  the 
only  corruption  to  the  SOI  is  receiver  noise.  The  third 
objective,  which  is  met  in  part  I,  is  to  introduce  multiple- 
sensor-pair  counterparts  of  some  of  these  various  TDOA- 
estimation  methods.  These  counterparts,  which  accom¬ 
modate  multiple  closely  spaced  sensors  as  well  as  widely 
spaced  platforms  and  both  wide-band  and  narrow-band 
signals,  suggest  that  the  dichotomy  between  the  array- 
based  direction  finding  problem  and  the  time-difference 
location  problem  is  not  as  strong  as  one  might  have 
thought. 

Since  the  new  methods  derive  their  superior  perfor¬ 
mance  by  exploiting  a  special  statistical  property  shared 
by  most  communications  and  telemetry  signals,  the  pre¬ 
sentation  begins  in  Section  II  by  describing  this  special 
property,  called  cyclostationarity,  in  a  suitable  theoretical 
framework.  This  leads  to  a  general  cyclostationary  model 
for  the  TDOA  problem  in  Section  III.  Then  in  Section  IV, 
each  of  a  number  of  methods  in  the  new  class  are  derived 
and  their  principles  of  operation  are  explained.  In  part  II, 
after  a  brief  introduction  in  Section  I,  specific  algorithms 
and  practical  considerations  for  digital  implementation  are 
presented  and  discussed  in  Section  II.  In  Section  Ill,  the 
ability  to  eliminate  the  resolution-limitation  problem  is 
demonstrated  qualitatively  using  the  results  of  simula¬ 
tions,  and  in  Section  IV,  the  performance  is  quantitatively 
assessed  by  evaluating  the  MSE  of  TDOA  estimates  ob¬ 
tained  from  Monte  Carlo  simulations  for  a  variety  of  en¬ 


vironments  and  performance  degrading  conditions  includ¬ 
ing  single,  multiple,  narrow-band,  and  wide-band  SNOI, 
weak  SOI,  SOI  and  SNOI  with  nearly  the  same  cycle 
frequencies,  and  error  in  the  algorithms’  knowledge  of  the 
cycle  frequency  of  the  SOI.  Conclusions  are  drawn  in 
Section  V. 

II.  Cyclostationarity  and  Spectral  Correlation 
Most  signals  encountered  in  communications  and  te¬ 
lemetry  systems  are  appropriately  modeled  as  cyclosta¬ 
tionary  time  series  rather  than  as  stationary  time  series. 
This  is  a  direct  result  of  underlying  periodicities  in  the 
time  series  due  to  periodic  sampling,  scanning,  modulat¬ 
ing,  multiplexing,  and  coding  operations  employed  in  the 
transmitter.  The  cyclostationarity  property  can  be  used  by 
receivers  to  obtain  substantial  tolerance  to  noise  and  in¬ 
terference  for  various  signal  processing  tasks  including 
TDOA  estimation  [10],  [13].  The  fundamentals  of  the 
spectral  correlation  theory  of  cyclostationary  time-series 
are  briefly  surveyed  in  this  section.  This  includes  the  limit 
cyclic  mean,  which  is  the  characteristic  property  of  what 
is  called  first-order  periodicity,  and  the  limit  cyclic  auto¬ 
correlation  and  limit  cyclic  spectrum,  which  are  the  char¬ 
acteristic  properties  of  what  is  called  second-order  peri¬ 
odicity.  The  spiectral  correlation  function,  an  alternative 
but  equivalent  characterization  of  second-order  periodic¬ 
ity,  is  also  described  and  its  equivalence  to  the  limit  cyclic 
spectrum  is  explained.  This  material  is  a  brief  summary 
of  the  comprehensive  treatment  given  in  [10]  that  is  pro¬ 
vided  for  the  reader’s  convenience.  Somewhat  more  de¬ 
tailed  surveys  are  given  in  [14],  [15],  [23]. 

A.  Cyclic  Autocorrelation  and  Cyclic  Spectrum 

Consider  a  time  series  x(u)  that  contains  an  additive 
periodic  component  with  period  To-  In  order  to  extract  the 
periodicity,  a  technique  called  synchronized  averaging  can 
be  used.  In  this  technique,  for  each  location  t,  2N  +  1 
points  are  sampled  from  x{u)  with  sampling  interval  Tq 
and  then  summed  to  obtain  the  time-variant  mean 

1  ^ 

^  ,  I  ^  +  ”^o)  (1) 

2N  +  \  n=-N 

where  T  =  {2N  -I-  1)  Tq  is  the  total  length  of  the  data  seg¬ 
ment  over  which  averaging  is  carried  out  for  each  value 
of  t.  As  the  number  of  samples  averaged  increases  without 
bound,  randomness  is  completely  averaged  out  and  the 
limit  time-variant  mean  is  obtained 

1  ^ 

Af,(r)  =  lim  Af,(r)r  =  lim  *  D  x{t  -I-  nTo). 

r-oo  JV-ao  2N  -r  1  /!=  -N 

It  can  be  shown  from  (2)  that  (t)  is  periodic  with  period 
To 

MAt  +  To)  =  MAt)  (3) 
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and  this  limit  function  is,  therefore,  called  the  limit  pe¬ 
riodic  mean.  Applying  the  Fourier  series  expansion  yields 
00 

M,{t)  =  S  (4a) 

m  ~  -  OD 

where  the  Fourier  coefficients  are  given  by 
I 

^^/To  A  i  (4b) 

M)  •-  -Tu/; 

Substituting  (2)  into  (4b)  produces  the  equivalence 
1 

=  Urn  -  .r(r)c"^”"'/^»dr 

r-oo  T  J-r/2 

4  <.x(r)f-''™'/^»>  (5) 

which  is  called  the  limit  cyclic  mean.  Further  substituting 
(5)  into  (4a)  yields  the  synchronized  averaging  identity 

1  ^ 

M,(t)  =  lim  -  -  Z  x(t  +  nTo) 

AI  — 00  (2/V  +  1)  n  =  -N 
00  «r/2 

=  S  lim  -  \  x(r -f- (6) 

m=  -00  T-(x  T  J-r/2 

For  the  situation  in  which  x(r)  contains  multiple  incom¬ 
mensurate  periodicities,  the  representation  (4),  (5)  can  be 
generalized  to 

MAt)  =  2  (7a) 

a 

where 

A/?  4  <x(r)e''^^“'>  (7b) 

with  a  representing  all  harmonics  (integer  multiples)  of 
the  fundamental  frequencies  of  additive  periodicities  con¬ 
tained  in  x(t). 

We  say  that  a  time  series  contains  first-order  periodicity 
with  frequency  a  if  and  only  if  A/“  is  not  identically  zero 
for  some  a  ^  0.  If  the  limit  cyclic  mean  Af“  exists  in  the 
temporal  mean-square  sense,  then  the  time  series  exhibits 
a  spectral  line  at  frequency  a  with  strength  |A/“p. 

For  a  time  series  x(,t)  that  does  not  contain  first-order 
periodicity  (spectral  lines),  but  does  contain  the  more  sub¬ 
tle  form  of  periodicity  called  second-order  periodicity,  or 
cyclostationarity,  synchronized  averaging  of  the  lag  prod¬ 
uct  can  be  used  to  extract  the  periodicity.  Consider  the 
time-variant  autocorrelation 

1  ^ 

2yV  -1-  1  n  =  -V 

■x*(t-  xjl  nTo)  (8) 

with  T  =  (2A1  -f-  l)ro.  Similar  to  (2),  the  randomness  is 
averaged  out  when  N  increases  without  bound,  and  the 
limit  time-variant  autocorrelation,  called  the  limit  peri¬ 
odic  autocorrelation,  is  obtained: 

1  '' 

RAt,  r)  =  lim  RAt,  t)t  =  Hm  —  ,  Z 

7* -*06  /V-»oo  2N  +  1 

•  x(t  •+  t/2  +  «7’o)JC*i,r  -  t/2  -I-  nTo). 


By  using  the  Fourier  series  expansion,  (9)  can  be  reex¬ 
pressed  as 

00 

RAt,  r)  =  Z  /?r''''“(T)e'‘”"^^‘'  (10a) 

m  =  '  OD 

where 

RUt)  =  <x(i  -I-  T/2)x*(r  -  T/2)e"‘^'“'>  (10b) 

which  is  called  the  limit  cyclic  autocorrelation,  and  the 
parameter  a  is  called  the  cycle  frequency.  For  a  =  0, 
(10b)  is  recognized  as  the  conventional  limit  autocorre¬ 
lation.  The  limit  cyclic  autocorrelation  is  a  characteristic 
property  of  second-order  periodicity,  or  cyclostationarity, 
in  that  x(t)  is  said  to  contain  second-order  periodicity  if 
and  only  if  the  limit  cyclic  autocorrelation  is  not  identi¬ 
cally  zero  for  some  nonzero  cycle  frequency  a.  When  a 
time-series  x(t)  contains  multiple  incommensurate  sec¬ 
ond-order  periodicities,  (10a)  generalizes  to 

/?,(/,  r)  4  Z/?“(r)^'^“'  (11) 

a 

with  a  representing  all  harmonics  of  the  fundamental  fre¬ 
quencies  of  second-order  periodicities  contained  in  x(t). 

The  limit  cyclic  spectrum  is  defined  to  be  the  Fourier 
transform  of  the  limit  cyclic  autocorrelation 

SCO 

/?,“(T)e-‘^’^dT.  (12) 

—  OD 

It  follows  from  (10b)  and  (11)  that  for  a  =  0  the  limit 
cyclic  spectrum  becomes  the  conventional  limit  spectrum 
(or  power  spectral  density). 

For  x(t)  real,  R“(t)  has  even  symmetry  in  r  and 
5“(/)  has  even  symmetry  in  /.  Also  /?“(r)  and  5“(/) 
have  conjugate  (Hermitian)  symmetry  in  a.  When  a  =  0, 
we  know  that  the  conventional  limit  autocorrelation  func¬ 
tion  has  the  property  f?®(0)  >  /?®(t).  This  property  gen¬ 
eralizes  to /?°(0)  >  R“(t). 

The  discrete-time  counterparts  of  the  limit  cyclic  au¬ 
tocorrelation  and  the  limit  cyclic  spectrum  are  given  by, 
respectively, 

1  ^ 

^^(iTA  ^  lim  — —  Z  x(nT,  +  iTAx*(nTA 

N-’xIN  +  1  B  =  -N 

■  exp  [-/2ira(n  +  k/2)TA  (13) 

and 

OD 

S"(f)  i  Z  (14) 

k  =  -00 

where  T,  denotes  the  sampling  period. 

By  applying  the  synchronized  average  identity  (6)  to 
the  lag  product  ofx(r),  the  relationship  b<*»ween  the  li’^it 
cyclic  autocorrelations  for  jc(/)  and  {jclnT,)}  is  found  to 
be 

00 

R^ikTA  =  Z  R°^'”^^AkTAe‘''^. 


(9) 


(15) 
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By  substituting  (15)  into  (14),  we  obtain  the  relationship 
between  the  two  spectral  correlation  functions; 


Finally,  letting  7  «  to  obtain  unlimited  spectral  reso¬ 

lution  results  in 


5?(/)  =  ^  £  -7).  (16) 

/,  n.m=-0D  \  27  TJ 

It  can  be  seen  that  5“  ( /)  is  a  periodic  function  of / and  a 
with  periods  of  1/7,  and  2/7^,  respectively;  that  is. 


1  f  f  r ,  1 

1  U™  + 


lim  lim  /)^ 

T-*  OB  il/  -•  00 


{/  +  7)  = 


In  addition,  5“(/)  is  jointly  periodic  in  /and  a;  that  is, 
(/- ^)  =  5?(/)  (17c) 

for  all  integers  n. 

B.  Spectral  and  Temporal  Coherence 

An  alternative  way  to  characterize  second-order  peri¬ 
odicity  is  in  terms  of  spectral  correlation.  Consider  the 
time-variant  finite-time  complex  spectrum  of  x(t) 

pf  +  r/2 

XT(t,f)^\  x{u)e-‘^^^du  (18) 

J/  -  r/2 

which  is  a  measurement  of  temporally  local  frequency 
content  of  xit)  in  the  interval  {/  -  7/2,  f  +  7/2]  with 
spectral  resolution  width  1/7.  The  complex  spectrum 
XAt,  f)  is  shifted  in  frequency  from  / to /  +  a/2  and  / 
-  a/2,  and  the  bandwidth-normalized  temporal  correla¬ 
tion  of  the  two  shifted  time-variant  complex  spectra  is 
measured, 

1  1 

,  -XtIsJ+oc/I) 
at  J(-^/2  7 

■  X^{s,f  -  a/2)  ds,  (19) 

Then  letting  At  -►  00  to  remove  ail  randomness,  (16)  be¬ 
comes  (using  (18)) 

lim 


-  U.  f'  ,  -  a 

r-t»  J-r/2  I  I 

j.r/2 

=  lim  1  Rj{T)e~‘'‘^  dr 

r-00  J-r/2 


-nrfr 


=  5;(/).  (21) 

Thus,  it  is  seen  from  (21)  that  the  limit  cycle  spectrum 
defined  by  (10b)  and  (11)  can  be  obtained  by  first  mea¬ 
suring  the  correlation  of  two  spectral  components  at  fre¬ 
quencies  f  +  a/2  and  f  -  a/2  and  then  idealizing  the 
measurement  by  letting  A/  -*  00  followed  by  letting  7  -* 
00.  Consequently,  the  limit  cyclic  spectrum  is  called  the 
spectral  correlation  function. 

The  spectral  correlation  function  can  be  reinterpreted 
as  the  conventional  cross  spectral  density  for  the  pair  of 
time  series 

u(0  =  x(/)e''’“'  (22a) 

v(t)  4  x(t)e‘’‘^  (22b) 

obtained  by  frequency  shifting  jr(0.  Since 

/?2r(r)  ^  <u(t  +  T/2)v*(t  -  t/2)>  =  /?“(T)  (23) 


5°.(/)^(  <,(T)e-'^^dr  =  S?(/).  (24) 

J  ~ao 

In  other  words,  the  density  of  correlation  between  spectral 
components  at  frequency  /  in  u(t)  and  v{t)  is  identical  to 
the  density  of  correlation  between  spectral  components  at 
frequencies /  -1-  a/2  and /  -  a/2  in  x(r).  Also,  since  it 
follows  from  (22)  that 

^«(/)  =  5?(/+ a/2)  (25a) 

5?,(/)  =  Sx(/- a/2)  (25b) 

then  the  normalized  cross-spectral  density 

A  _ ^uv(f) 


c„xf)  4 


tS“(/)S?,(/)]' 


which  is  the  spectral  cross-coherence  function,  is  identi¬ 
cal  to 

n  I  t\  —  _ ^4  if)  A  r\ 


"MS  *"*“  "  dudv 

-r/2 

=  r  /?“(r)fl  - 
J-r/2  L  T] 


which  is  called  the  spectral  autocoherence  function  [10]. 

The  spectral  autocoherence  function  is  a  frequency  de¬ 
composed  correlation  coefficient  (assuming  that  = 
0  so  that  the  correlation  5“(/)  is  a  covariance)  for  the 
frequency-shifted  versions  (22)  of  Jt(r).  Thus,  it  has  the 


GARDNHR  AND  CHKN  SIGN  AL  SbLbX  I  IVK  FDOA  HSTIMATIUN,  PART  I 


1173 


property  |C‘‘(/)|  ^  1,  which  is  held  by  all  correlation 
coefficients.  Analogously,  the  quantity 


KAt) 


(28) 


which  is  called  the  temporal  autocoherence  function,  is 
the  correlation  coefficient  (assuming  that  m'*  =  0)  for  the 
time-shifted  versions 


u'(t)  ^  xU  +  tID  (29a) 

v'(t)  =  .r(/  -  t/2)  (29b) 


of.t(/).  Similarly,  the  more  general  temporal/spectral  au¬ 
tocoherence  function' 


R'I(t) 

/??(0) 


is  the  correlation  coefficient  (assuming  /m“ 
time-  and  frequency-shifted  versions 

u"(t)  =  x(t  ^  t/2) 

v'At)  =  x(t  -  rlDe 


(30) 
0)  for  the 

Ola) 

(31b) 


of  jc(r).  The  Fourier  transform  of  the  unnormalized  ver¬ 
sion  /?"(t)  of  /f"(r)  is  the  unnormalized  spectral  autoco¬ 
herence  function  S"(/). 

For  a  pair  of  time-series  x{t)  and  y(r),  the  natural  gen¬ 
eralization  of  /f“(T)  is  given  by 


r;(T)  = 


RUr) 


(/??(0)«?(0)1'^- 

and  the  natural  generalization  of  C°(f)  is  given  by 

C.,n= _ _ 

|S;(/+»/2)S?(/-c/2)|'” 

where 


(32) 


(33) 


RUt)  ^  <yU  +  T/2)x*(t  -  r/2)e-'^"“'>  (34) 


and 


Sl(f) 


(35) 


or 

S"Af)  =  lim  lim 

T  -*  00  A!  -*  oi 

in  which 


At 


Jl  -t-Alf  2 
l-Al/2 


-Yris,f+  a/2) 


(36) 


■  XUsJ  -  a/2)  ds.  (37) 

Conventional  methods  of  TDOA  estimation,  some  of 
which  are  described,  for  example,  in  [2],  [11],  and  [16], 


"Although  it  would  seem  to  be  consistent  to  call  R"{0)/R°(0)  the  spectral 
autocoherence  function,  this  term  is  reserved  for  (27). 


are  based  on  the  unnormalized  temporal  cross  coherence 
function  R^Ax)  and  its  Fourier  transform  S^’,(  /).  The  new 
methods  of  TDOA  estimation  presented  in  Section  IV  aic 
based  on  the  unnormalized  temporal/spectral  cross  coher¬ 
ence  function  R‘'Ax)  and  its  Fourier  transform  S“A  /)■ 
unnormalized  spectral  cross  coherence  function.  Exam¬ 
ples  of  the  functions  R\‘Ax)  and  S/A  f)  (o'"  numerous  types 
of  communication/telemetry  signals  are  described  in  [  10], 
[13],  [17],  ]18].  One  such  example  is  given  in  the  Ap¬ 
pendix.  As  explained  in  the  references,  typical  cycle  fre¬ 
quencies  a  in  these  signals  include  the  keying  rate,  a 
cr*,  and  its  harmonics  in  keyed  digital  systems,  such  as 
amplitude-,  phase-,  and  frequency-shifted  keying  (ASK, 
PSK,  and  FSK),  and  the  doubled  carrier  frequency,  a  = 
2/,  ,  in  continuous-wave  analog  systems,  such  as  ampli¬ 
tude,  phase,  and  frequency  modulation  (AM,  PM,  and 
FM),  as  well  as  sums  and  differences  of  these  cycle  fre¬ 
quencies  in  digital  carrier-modulated  systems  such  as 
ASK,  PSK,  and  FSK. 

Because  of  the  current  and  growing  popularity  of  PSK 
in  both  communications  and  telemetry  systems,  the  sim¬ 
ulations  in  Part  II  of  this  paper  focus  on  this  type  of  signal 
for  the  SOI,  but  consider  both  PSK  and  AM  for  the  SNOI. 
Since  carrier  frequencies  used  in  the  transmitter  are  easily 
changed,  they  can  be  more  difficult  for  an  unintended  re¬ 
ceiver  to  have  knowledge  of  them.  Consequently,  only 
the  keying  rate  in  the  PSK  SOI  is  used  in  the  simulations 
presented  in  Part  11.  Nevertheless,  it  is  pointed  out  that 
for  signals  that  do  possess  the  cycle  frequency  a  =  2/ , 
the  spectral  correlation  at  a  =  2/.  is  typically  stronger 
than  that  at  a  =  a*  and,  consequently,  the  performance 
of  TDOA  estimators  that  exploit  a  =  2/  is  typically  su¬ 
perior.  It  also  should  be  clarified  that  when  complex  data 
x(t)  is  used,  and  a  equal  to  or  associated  with  a  carrier 
frequency  is  used,  then  the  conjugate  cyclic  correlations, 
for  which  the  conjugate  in  (10b)  and  (34)  is  removed,  must 
be  used.  As  a  consequence  of  this,  the  conjugates  in  (19) 
and  (37)  must  be  removed  from  the  second  factor  and  the 
argument  /  -  a/2  in  this  factor  must  be  changed  to  a/2 
-/[lO]. 


III.  Cyclostationary  Model  for  TDOA 

Before  proceeding  to  the  introduction  of  the  cyclosta- 
tionarity-exploiting  methods  for  estimating  TDOA,  let  us 
first  determine  the  general  form  of  the  functions  /?“  (r) 
and  5“  (/),  as  well  as  /?“(t).  R/(t),  S/if),  and  S/Af), 
for  the  particular  model  of  xit)  and  y(r)  that  describes  the 
data  received  from  a  pair  of  antenna  elements. 

A  signal  (SOI)  radiating  from  a  remote  source,  propa¬ 
gating  through  a  nondispersive  medium,  and  impinging 
on  two  spatially  separated  antenna  elements,  together  with 
interfering  signals  (SNOI),  and  received  in  the  presence 
of  receiver  (including  antenna)  noise,  can  be  modeled  as 

x{t)  =  5(r)  +  nit)  (38a) 

y(r)  =  Asit  -  D)  +  mit)  (38b) 


I  I  u 
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where  s(t)  is  the  SOI,  n{t)  and  mU)  are  the  SNOI,  includ¬ 
ing  noise,  D  is  the  TDOA  to  be  estimated  using  measure¬ 
ments  x(f)  and  >’(/),  and  A  represents  magnitude  (and  phase 
for  complex  data)  mismatch  between  receivers.  In  this 
ideal  model  all  signal  distortion  due  to  the  propagation 
channels  and  the  sensors  is  assumed  to  be  negligible  or  to 
be  matched  (i.e.,  matched  sensors  and  matched  chan¬ 
nels).  The  effects  of  mismatch  here  are  described  in  part 
II  (!].  It  is  also  assumed  that  s{t),  n(t),  and  m{t)  have  zero- 
mean  (time  average)  values,  and  that  s(t)  is  statistically 
independent  (over  time)  of  n(t)  and  m(t).  However,  since 
n{t)  and  m(t)  can  contain  the  same  interfering  signals  (with 
different  times  of  arrival),  as  well  as  independent  receiver 
noise,  then  they  can  be  statistically  dependent.  It  follows 
that  the  limit  cyclic  cross  correlation  and  autocorrelations 
are  given  by 


/?“  (T)  =  AR^T  -  Die-""'’  +  /?L(r) 

(39a) 

/?“(T)  =  /??(T)  +  RUt) 

(39b) 

/?“(r)  =  |/l|'/?r(r)e +  /?“(t). 

(39c) 

The  limit  cyclic  cross  spectrum  and  autospectra  can  be 
obtained  by  Fourier  transforming  (39),  and  are  given  by 

S“(/)  =  ^S?(/)e-'2'< Sl(/) 

(40a) 

SUf)  =  5?(/)  +  S:(f) 

(40b) 

S“(/)  =  l4pS“(/)e S“(/). 

(40c) 

Although  these  ideal  (infinite-time  average)  measure¬ 
ments  contain  the  TDOA  parameter  D  to  be  estimated,  the 
parts  of  these  measurements  that  contain  the  parameter  D, 
namely,  those  dependent  on  the  SOI  s(t),  are  masked  by 
the  parts  due  to  the  SNOI  n(t)  and  m(t).  As  a  result,  when 
the  SNOI  are  both  temporally  and  spectrally  coincident 
with  the  SOI  (and  knowledge  required  to  do  spatial  filter¬ 
ing  to  reject  the  interfering  signals  in  n(r)  and  m(t)  prior 
to  recording  jc(r)  and  y(r)  is  not  available),  then  conven¬ 
tional  methods  which  use  these  measurements  with  only 
a  =  0  cannot  avoid  the  resolution-limitation  problem  de¬ 
scribed  in  Section  I.  Furthermore,  even  if  the  individual 
cross-correlation  peaks  due  to  the  SOI  and  SNOI  can  be 
resolved,  the  conventional  methods  still  have  to  deter¬ 
mine  which  of  the  multiple  peaks  is  due  to  the  SOI  and 
which  are  due  to  the  SNOI. 

In  contrast  to  this  difficult  situation,  if  the  SOI  exhibits 
a  cycle  frequency  a  not  shared  by  any  of  the  SNOI  (this 
is  usually  the  case  since  keying  rates  and  carrier  frequen¬ 
cies  of  a  SOI  are  rarely  the  same  as  those  of  the  SNOI 
even  though  signals  reside  in  the  same  band),  then  by 
using  this  value  of  a  in  the  measurements  (39),  (40),  we 
obtain  (ideally,  for  infinite  averaging  time) 

/?“(T)  =  /?“(r)  s  /?“„(T)  s  0  (41) 

and,  therefore, 

5“(/)  =  SZ(f)  ^  ^  0.  (42) 


This  is  the  means  by  which  we  obtain  signal  selectivity  in 
our  measurements.  Using  (41)  and  (42),  (39)  and  (40)  re¬ 
duce  to 


Rlir) 

=  A/?“(t  -  D)e~'^^ 

(43a) 

R^t) 

=  Rlir) 

(43b) 

R^t) 

=  |/lp/?f(T)e-'‘'“" 

(43c) 

SUf) 

(44a) 

S“(/) 

=  5“(/) 

(44b) 

Sv(/) 

=  iApS“(/)e-'''“'’. 

(44c) 

These  ideal  measurements  are  the  same  as  those  we  would 
obtain  directly  from  the  idealized  SNOI-free  model 

x(t)  =  s(0  (45a) 

y(f)  =  As(t  -  D).  (45b) 

Observe  that  these  same  results  hold  true  when  there 
are  multiple  SOI  in  the  environment,  each  exhibiting  a 
distinct  cycle  frequency  a.  In  practice,  if  there  are  signals 
of  potential  interest  whose  cycle  frequencies  a  are  not 
known  in  advance,  these  cycle  frequencies  can  easily  be 
estimated  prior  to  the  execution  of  a  TDOA-estimation 
algorithm  that  requires  knowledge  of  a.  Specifically,  all 
cycle  frequencies  of  all  received  signals  can  be  estimated 
by  simply  seeking  spectral  lines  in  the  FFT’s  of  one  or 
more  lag-product  waveforms  z,(r)  =  x(t  -  T)x*{t)  cor¬ 
responding  to  one  or  more  lag  values  r  [1].  (Also,  one 
can  do  a  least  squares  fit  over  t  of  z,(0  to  for  all  a 
of  potential  interest.)  This  capability  is  demonstrated  in 
[19],  where  it  is  used  in  conjunction  with  high-resolution 
array-based  algorithms  that  exploit  cyclostationarity  to 
obtain  signal  selectivity  for  direction  finding.  Also,  the 
sensitivity  of  the  performance  of  the  TDOA-estimation 
methods  introduced  in  this  paper  to  errors  in  a  is  quanti¬ 
tatively  evaluated  in  part  II  [  I  ] . 

In  the  next  section,  we  use  these  SNOI-free  relations 
(43)-(45)  to  design  a  number  of  algorithms  for  estimating 
the  TDOA  with  tolerance  to  SNOI. 

IV.  Cyclostationarity-Exploiting  Methods  tor 
TDOA  Estimation 

A.  Conventional  GCC  Methods 

The  conventional  cross-correlation  method  for  TDOA 
estimation  is  based  on  (43a)  with  a  =  0,  which  applies 
only  in  the  absence  of  SNOI.  Let  us  denote  this  idealized 
measurement  by 

^  R%(r)  =  /l/f?(T  -  D).  (46) 

Since  the  limit  autocorrelation  /?”(«)  peaks  at  «  =  0  for 
all  signals,  then  peaks  at  the  TDOA  t  =  D.  Thus, 
the  value  of  r  at  which  an  estimate  (obtained  from  finite 
time  averaging)  of  the  ideal  OoIt)  peaks  can  be  taken  as 
an  estimate  of  D. 

Conventional  generalized  cross-correlation  methods  for 
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TDOA  estimation  can  be  based  on  (44a)-(44c)  with  a  = 
0,  which  applies  only  in  the  absence  of  SNOI: 


SUf)  =  AS'^Af)e-‘~^^‘^ 

(47a) 

5?(/)  =  5?(/) 

(47b) 

SUf)  =  |/l|-5?(/). 

(47c) 

For  example,  if  the  support  band  of  S°(/)  is  as  follows: 

'>0,  ll/l  -/ol  <  Bo/2 

=  0,  otherwise 


sUf)  = 


then  the  ratio  of  (ideal)  measured  spectra  is  simply 
SUf)  ^(Ae  1 1  /I  -  /„  I  <  Bo/2 

S^x(f)  (o,  otherwise 


(where  the  value  zero  is  assigned  by  the  TDOA  estimation 
method  since  there  is  no  relevant  information  in  this  ratio 
outside  the  support  band  and  0/0  is  undefined).  Inverse 
Fourier  transforming  this  ratio  of  spectra  yields  the  ideal¬ 
ized  measurement 


bo(.T) 


sUf) 

I/I -/,!<&/’  S°(f) 


iTt/t 


df 


(48a) 


A  sin  [irBo(T  —  D)} 
7r(T  -  D)/2 


cos  [2ir/o(r  -  D)] 


(48b) 


which  peaks  at  t  =  D,  as  desired.  Because  of  the  weight¬ 
ing  of  the  cross-spectrum  5?,(/)  by 


mf) 


r  1 
.  s°{fr 
io, 


ll/I  -/ol  <  Bo/2 
otherwise 


the  TDOA-estimation  method  based  on  an  estimate  of 
bo(T)  is  called  a  generalized  cross-correlation  (GCC) 
method.  With  W(f)  =  1,  it  reduces  to  the  cross-corre¬ 
lation  method  based  on  an  estimate  of  ao(T).  Other  GCC 
methods  simply  use  other  weighting  functions  H^(  /  ),  typ¬ 
ically  chosen  to  minimize  the  effects  of  noise  and  SNOI 
as  revealed  in  (39),  (40)  [2],  [11],  [16].  However,  since 
the  use  of  any  weighting  function  fV(f)  can  always  be 
reinterpreted  as  simply  applying  the  crosscorrelation 
method  based  on  an  estimate  of  ao(T)  to  spectrally  filtered 
data  j:(r)  and  >(/),  then  whenever  some  SNOI  eover  the 
same  spectral  band  as  the  SOI,  there  is  little  that  GCC 
methods  can  do  to  reduce  the  undesirable  effects  of  those 
SNOI. 


B.  CCC  Method 

To  circumvent  this  interference  problem  with  the  GCC 
method,  the  cyclic  counterpart  of  the  conventional  method 
based  on  an  estimate  of  ao(r)  in  (46),  which  is  called  the 
cyclic  cross-correlation  (CCC)  method,^  uses  an  estimate 


'The  CCC  method  was  invented  by  Gardner  and  first  disclosed  in  the 
1985  manuscript  for  the  book  (10),  which  appeared  in  1987.  However,  it 
was  also  proposed  in  December  of  1987  by  Friedman  et  al.  (20),  who  re¬ 
ceived  a  patent  for  an  implementation  of  it  in  December  1989. 


of  the  idealized  measurement  from  (43a) 

ttjr)  =  |B:,(r)|  =  MIlBrir  -  D)|  (49) 

which  is  valid  even  in  the  presence  of  SNOI.  Although 
the  peak  of  |B°(u)|  does  not  occur  at  u  =  0  for  all  signals, 
since  !B“(u)|  is  an  even  function  of  u,  there  is  typically  a 
pair  of  peaks  straddling  u  =  0  as  demonstrated  in  Fig. 
1(a)  for  the  binary  PSK  (BPSK)  signal  with  rectangular 
keying  envelope,  cycle  frequency  a  equal  to  the  keying 
rate  a^,  and  TDOA  D  =  487"^,  as  described  in  the  Appen¬ 
dix.  This  pair  of  peaks  appears  in  the  envelope  of  a  func¬ 
tion  that  oscillates  with  frequency  equal  to  the  carrier  fre¬ 
quency  of  the  signal.  Thus,  the  midpoint  of  the  pair  of 
peaks  in  the  estimate  d^ir)  of  a^ir)  can  be  taken  as  an 
interference-tolerant  estimate  of  the  TDOA  D. 

As  an  alternative,  the  estimated  CCC  function  /?“j(t) 
can  be  correlated  with  its  idealized  version  evaluated  at  D 
=  0,  namely,  /?“ (T),and  a  maximum  can  then  be  sought; 

D  4  arg  max  {d;(T)}  (50) 


where 


dL(T)  ^  |/?“  (r)  @  B?(t)| 


j  /?“  (T  +  «)/?“(«)♦  du 


(51) 


where  @  denotes  correlation.  We  refer  to  this  alternative 
(50)-(51)  as  the  correlated  CCC  (CCCC)  method.  By 
substituting  (43a)  into  (51),  we  obtain  the  ideal  CCCC 
function 


a;(T)  =  \A\  |B“(r  +  D)  ®  R°{t)\  .  (52) 

A  graph  of  this  ideal  CCCC  function  for  the  BPSK  signal 
specified  in  the  Appendix  with  rectangular  keying  enve¬ 
lope,  a  =  a^,  and  D  =  487’j  is  shown  in  Fig.  1(b),  where 
it  can  be  seen  that  there  is  a  unique  peak  (in  the  envelope) 
at  T  =  D; 


arg  max  {a;(r)}  =  D.  (53) 

T 

If  the  ideal  cyclic  autocorrelation  /?“(t)  used  in  (51)  is 
not  known  (to  within  a  complex  factor  that  is  independent 
of  t),  then  it  can  be  replaced  with  the  estimate  B“(r)  since 
B?(t)  =  B“(t)  (cf.  (43b)): 

d"(T)  4  1«“(t)@R“(t)|.  (54) 

The  ideal  version  o"  (t)  of  this  CCCC  function  is  exactly 
the  same  (within  a  scale  factor)  as  (52)  and,  therefore, 
peaks  at  the  TDOA  r  =  D.  This  CCCC  method'* 

D  =  arg  max  {d"(T)}  (55) 

T 

is  equivalent  to  the  solution  to  a  least  squares  fitting  prob¬ 
lem  as  explained  following  the  next  method. 

*The  CCCC  method  (54)  is  equivalent  to  the  SPECCOA  method  which 
was  invented  by  Gardner  in  1987. 
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Fig.  1 ,  (a)  The  ideal  CCC  function  a„(T)  in  (49)  with  o  =  a,  for  the  BPSK 
signal  with  rectangular  keying  envelope  and  TDOA  D  =  487',.  (b)  The 
ideal  CCCC  function  a'„{r)  in  (52)  with  a  =  a,  for  the  BPSK  signal  with 
rectangular  keying  envelope  and  TDOA  D  =  48T,. 


C  SPECCORR  Method 

The  CCC  method  can  be  substantially  improved  on 
through  the  use  of  generalized  CCC  methods.  For  exam¬ 
ple,  the  cyclic  counterpart  of  the  GCC  method  based  on 
an  estimate  of  boir)  in  (48a)  uses  an  estimate  b„(T)  of  the 
idealized  measurement 


.  ,  ,  A  [  SyAf)  -2,/^ 

^a(T)  —  \  aa/  ^ 


=  s, 


=  \A\ 


\f\-fa\<Bj2 

sin  1‘kBJt  -  P)1 
7r(T  -  D)/2 


(/) 

^^Cv/tr-D) 


(56a) 

(56b) 


cos  [2T/a(7  -  D)] 


(56c) 


where  /„  and  Ba  are  the  center  and  width  of  the  support 
band  for  5“(/)  (e.g.,  for  BPSK,  if  we  choose  a  =  2/^, 


we  have /„  =  0  and  =  By;  but  if  we  choose  a  =  a*, 
we  have /„  =  /.  and  B^  =  2at,  see  the  Appendix).  Like 
the  CCCC  method,  this  method 

D  arg  max  (57) 

T 

called  the  spectral  correlation  latio  (SPECCORR) 

method,^  has  the  advantage  over  the  CCC  method  that  it 

always  produces  a  peak  at  the  TDOA  D  (when  idealized— 
infinite-time  average— measurements  are  used)  19].  [lO); 

arg  max  {fia(r)}  =  D.  (58) 


When  knowledge  of  the  parameters  /„  and  is  avail¬ 
able,  we  call  this  the  band-limited  SPECCORR  (BL- 
SPF.CCORR)  method.  On  the  other  hand,  when/^  and  B^ 
are  not  available,  and  integration  in  (52a)  is  carried  out 
over  the  entire  reception  band  of  the  source  location  sys¬ 
tem,  then  the  modifier  BL  is  not  used.  When  and  are 
not  knov'n  in  advance,  they  can  be  selected  to  match  the 
band  over  which  the  magnitude  of  the  estimate  of  the  ratio 
in  the  integrand  in  (54a)  is  sufficiently  close  to  unity. 
Moreover,  within  this  band,  the  actual  magnitude  can  be 
replaced  with  unity  to  reduce  undesirable  random  effects. 

The  SPECCORR  method  can  be  arrived  at  through  a 
least  squares  optimization  procedure.  Specifically,  since 
the  ideal  spectral  correlation  ratio  is  given  by 

■5«(/) 

=  C  exp  ( -i(2ir(  f  +  a/2)D  -  «])  (59) 


where  C  =  1/4)  and  <f>  =  arg  {/4},  then  we  can  seek  a  least 
squares  fit  of  to  an  estimate  of  the 

spectral  correlation  ratio: 


min 

C. 


1 1 /I  -h\<B„n 


suf) 

5“(/) 


-  C 


■  exp  (-i[27r(/+  a/2)T  -  0]) 


(60) 


The  solution  is  given  by 

D  =  arg  max  {fi„(T)} 

r 

where  6a  (r)  is  an  estimate  of 


6a  (t) 


4 

V  I 


ll/|-/ai  -Ba/. 


5«(/) 


^“(/) 


,  i2  wfr 


df 


(61) 


(62) 


which  is  identical  to  (56a).  A  graph  of  the  ideal  BL- 
SPECCORR  function  (56c)  for  the  BPSK  signal  described 
in  the  Appendix  with  a  =  a*  and  D  =  487,  is  shown  in 
Fig.  2. 

Before  proceeding  to  discuss  other  generalized  CCC 
methods,  it  is  pointed  out  that  unlike  the  GCC  methods, 
which  just  use  real- valued  weighting  functions  W(f)  to 
emphasize  certain  spectral  bands  and  deemphasize  others. 


’The  SPECCORR  method  was  invented  by  Gardner  in  1986.  was  devel¬ 
oped  and  reduced  to  practice  jointly  by  Gardner  and  Chen,  and  first  ap¬ 
peared  in  (10]  in  1987. 
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Fig.  2,  (a)  The  ideal  leasi  squares  BL-SPECCORR  funclion  b,.(r)  in  (56c) 
for  complex  data  with  a  =  a,.  B„  =  la,,,  and/,  =  /  for  the  BPSK  signal 
with  TDOA  D  =  dST,.  (b)  The  ideal  least  squares  BL-SPECCORR  func¬ 
tion  for  real  data  with  a  =  a^.  B„  =  2aj,  and/,  =  /  for  the  BPSK  signal 
with  TDOA  D  =  47r.. 


the  weighting  functions  in  the  generalized  CCC  methods 
are  complex  valued  and  contain  phase  information  about 
the  SOI.  Thus,  they  do  much  more  than  just  spectral 
weighting  in  the  usual  sense. 


where  c^lr)  is  an  estimate  of 

fa(r)  ^  J  |/?“(U)  -  C/??(U  -  r)c I' du 

=  J  -  r)|' 


The  solution  is  given  by 

D  =  arg  max  {c;(t)} 

r 

where  c4(t)  is  an  estimate  of 

c'Jt)  =  I  /?“(u)/?“{u  -  t)*  du 

SUf)SUf)*e‘-^^  df 


-\ 


(65) 

(66) 


(67a) 

(67b) 


and  where  (67b)  follows  from  (67a)  by  using  Parseval’s 
relation  for  the  Fourier  transform.  This  method  can  be 
interpreted  in  terms  of  maximizing  the  correlation  over / 
ol  two  spectral  correlation  functions  by  phase  alignment 
through  adjustment  of  the  linear  phase-versus-frequency 
term  2ir(/  +  ci/2)t.  Consequently,  this  approach  is 
called  the  spectral  coherence  alignment  (SPECCOA) 
method’  [12].  By  analogy  with  the  SPECCORR  method, 
it  is  also  called  the  spectral  correlation  product  (SPEC- 
CcRP)  method  to  distinguish  it  from  SPECCORR,  which 
also  has  the  SPECCOA  interpretation  1 13]. 

To  see  that  c^(t)  does  indeed  peak  at  t  =  D.  we  simply 


substitute  (44a).  (44b)  into  (67b)  to  obtain 

«r)  -  C 

J  |S?(/)|’exp  [/2t(/+  a/2)(T 

-  Z?)]  df 

(68a) 

» 

\S^(f)\'df 

(68b) 

=  c;(D). 

(68c) 

D.  SPECCOA  Method 

Other  generalized  CCC  methods  can  be  obtained  using 
various  ad  hoc  least  squares  optimization  procedures.*  For 
example,  since  (43a)  and  (43b)  reveals  that  (ideally) 

/?“(«)  =  CR^(u  -  D)e (6 

then  we  can  seek  the  value  of  the  estimate  t  =  D  that 
minimizes  the  sum  (over  the  lag  parameter  «)  of  squares 
of  error  magnitudes  between  the  measurements  of  the  left 
and  right  sides  of  (63)  with  D  replaced  by  t: 

D  =  arg  min  {4(t)}  (64) 


Therefore, 

arg  max  {c;(t)}  =  D.  (69) 

T 

A  graph  of  the  ideal  SPECCOA  function  c^ir)  given  by 
(68a)  is  shown  in  Fig.  3  for  the  BPSK  signal  specified  in 
llic  Appendix  with  a  =  a*  and  D  =  48T,. 

An  advantage  of  the  SPECCOA  function  c^It)  over  the 
SPECCORR  function  b„  (t)  is  that  the  former  contains  the 
factor  iS“(/)|’  in  the  integrand  (cf.  (68a)),  whereas  the 
latter  does  not  (cf  (56b)).  This  factor  can  deemphasize 
the  integrand  in  spectral  bands  where  there  is  little  con¬ 
tribution  from  the  SOI,  without  prior  knowledge  of  these 


^Although  clas.xical  statistical  principles,  such  as  maximum  likelihood, 
provide  an  alternative  approach,  we  have  not  yet  found  such  approaches  to 
be  tractable  for  the  non-Gaussian  nonstationary  models  of  interest. 


’The  SPECCOA  method  was  invented  by  Gardner  in  1987,  was  devel¬ 
oped  and  reduced  to  practice  jointly  by  Gardner  and  Chen,  and  was  first 
publicly  disclosed  at  a  1988  workshop  (12). 
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(b) 


Fig.  3.  (a)  The  ideal  lea.sl  squares  SPECCOA  function  c.'.ir)  in  (69a)  for 
complex  data  with  «  =  Oi  for  the  BPSK  signal  with  half-cosine  keying 
envelope  and  TDOA  D  —  487',.  (b)  The  ideal  least  squares  SPECCOA 
function  for  real  data  with  a  =  a*  for  the  BPSK  signal  with  half-cosine 
keying  envelope  and  TDOA  D  =  487,. 

bands.  As  the  simulations  in  part  II  show,  SPECCOA  can 
indeed  require  less  data  for  a  MSE  of  the  estimate  D  that 
is  comparable  to  that  of  SPECCORR  and  even  BL-SPEC- 
CORR! 

It  is  noted  that  SPECCOA,  as  expressed  in  (67a)  (after 
using  the  change  of  variables  u  -  t  =  v),  is  identical  to 
the  ad  hoc  method  CCCC  in  (54).  When  real-valued  data 
is  used,  the  scalar  A  in  the  model  (38)  is  real  and,  as  a 
result,  the  solutions  to  the  least  squares  problems  (60)  and 
(64)  with  v’  =  0  use  the  real-part  operation  Re  { • }  in 
place  of  the  magnitude  operation  |  •  |  in  (62)  and  (67).  The 
corresponding  SPECCORR  and  SPECCOA  functions  for 
real  data  are  shown  in  Figs.  2(b)  and  3(b),  respectively. 

E.  SPECCON  Method 

Another  ad  hoc  optimization  procedure  is  to  minimize 
with  respect  to  r  the  sum  (over  the  lag  parameter  u)  of 
magnitude  squared  cyclic  autocorrelation  values  for  the 
difference  signal  z(t)  =  Ax(t)  -  y(t  -F  t): 

/??(«)  =  |/1T/?“(m)  + 

-  ARUu  -  7)^'"“^  -  /!*/?“(«  +  T)e'’^“'  (70) 


where,  for  complex  data,  A  is  complex.  This  is  motivated 
by  the  fact  that  this  cyclic  autocorrelation  would  ideal.y 
be  identically  zero  if  and  only  if  t  =  D  and  A  =  A,  as 
revealed  by  substitution  of  (43a)-(43c)  into  (70).  Since 
Parseval’s  relation  yields 

J  |R“(«)|'’  du  =  J  \S‘^{f)\- df  =  djT)  (71) 

then  minimization  of  this  sum  of  magnitude-squared  val¬ 
ues  yields  the  spectral  coherence  nulling  (SPECCON) 
method* 

D  -  arg  min  {da(T))  (72) 

A.  7 

where  d„(T)  is  an  estimate  of  d^iT). 

For  applications  in  which  A  =  1  is  an  adequate  esti¬ 
mate,  by  substituting  (70)  into  (71)  and  eliminating  terms 
that  are  independent  of  t  we  obtain  the  following  equiv¬ 
alent  TDOA  estimate: 

D  =  arg  min  {d'„(T)}  (73) 

T 

where  d'^ir)  is  an  estimate  of 

J  [5“(/)5?(/)*f> 

+  S"(/)*S“.(/)e-'"'^^ 

-  [5?(/)*S“,(/)  -F  S"(/)S“(/)*) 

.  g  -i2t(  f- a/2)T 

-  [5“(/)5“(/)*  +  S“(/)*5“  (/)] 

.  ^-,2»(/+a/2)rj  (74) 

By  substituting  (44a)-(44c)  with  A  =  1  into  (74),  we  ob¬ 
tain  the  ideal  SPECCON  function 

d'jT)  =  Re  J  |5?(/)|^  [^--2T0,(r-D) 

_  2g  '2i(/-c»/2)(t-D)  _  2^ -i2»(/+a/2)lT-/»  J  ^ 

(75) 

A  graph  of  this  function  for  the  BPSK  signal  specified  in 
the  Appendix  with  a  =  a*  and  D  =  487’^  is  shown  in  Fig. 
4(a).  For  this  BPSK  signal,  it  can  be  seen  that  d'„(T)  not 
only  exhibits  a  global  minimum  at  t  =  D,  but  also  ex¬ 
hibits  local  minima  on  both  sides  of  D.  For  finite-time 
averaging,  the  actual  global  minimum  can  occur  at  one  of 
the  locations  of  the  ideal  local  minima.  Because  these  lo¬ 
cal  minima  are  symmetric  about  the  ideal  global  mini¬ 
mum,  a  pattern  recognition  technique  can  overcome  the 
problem  of  choosing  the  wrong  minimal  point  for  the 

"The  SPECCON  method  was  invented  by  Gaidner  in  1988.  was  devel¬ 
oped  and  reduced  to  practice  jointly  by  Gardner  and  Chen,  and  was  first 
disclosed  by  Chen  in  [21]. 
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Lig.  4.  (al  The  ideal  SPECCON  t'unelinn  </,',(r)  in  115)  with  «  =  ai  lor 
ihe  BPSK  Mgnal  wilh  hall-eosine  keving  envelope  and  TDOA  D  =  48r,. 
(h)  The  ideal  PP  SPECCON  lunelion  </"(rl  eorresponding  (o  (761  (using 
(7S|  in  plaee  of  i/  .iTlI  wilh  «  =  n,  for  Ihe  BPSK  signal  with  half-eosine 
keying  envelope  and  TDO.A  D  =  48/’,, 


TDOA  estimate  D.  Another  approach  is  to  circumvent  the 
multiple  minima  by  correlating  the  estimated  SPECCON 
function  with  the  ideal  function  d'„{r)  evaluated  at 
D  =  0  and  truncated  to  the  interval  |t|  <  1  /Za*  (to  re¬ 
move  excessive  oscillation),  and  then  seeking  a  maximum 
with  respect  to  t.  A  graph  of  the  ideal  version  (with 
d!,{T)  used  in  place  of  d’,(T))  of  this  postprocessed  SPEC¬ 
CON  (PP-SPECCON)  function 

d"(T)  =  d'„(T)  ‘i!  d'„(T)  (76a) 

where 

|r/,',(T)lD-ii  kl  ^  1/20!* 

0.  1t|  >  \/2a,  (76b) 

is  shown  in  Fig.  4(b)  for  the  BPSK  signal  with  a  =  a* 
and  D  =  487,.  It  can  be  seen  that  it  exhibits  a  unique  peak 
at  T  =  D.  as  desired. 

Although  the  simulations  presented  in  part  II  for  A  = 
I  show  that  the  SPECCON  method  (73).  without  pattern 
recognition,  performs  comparably  to  the  BL-SPECCORR 
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method,  they  also  show  that  the  PP-SPECCON  method 


D  =  arg  max  {^"(t)}  (77) 


performs  substantially  better  than  even  the  SPECCOA 
method  in  that  it  yields  comparable  MSE  with  less  data. 
Nevertheless,  it  does  require  knowledge  of  the  ideal  spec¬ 
tral  correlation  magnitude  function  |S‘'(/)|  for  the  SOI. 
whereas  SPECCOA  does  not  require  any  knowledge  about 
the  SOI,  except  of  course  a  cycle  frequency  a.  Also,  the 
amount  of  computation  required  by  SPECCON  (even 
without  postprocessing)  is  about  six  times  that  required 
by  SPECCOA.  Moreover,  this  simplified  version  of 
SPECCON  and  the  conesponding  PP-SPECCON  are  not 
appropriate  when  A  in  (38b)  is  not  close  to  unity. 

Although  one  might  expect  a  postprocessed  version  of 
SPECCOA,  analogous  to  (76),  to  outperform  SPECCOA, 
this  has  not  been  found  to  be  the  case  in  any  of  the  sim¬ 
ulations  performed  for  a  BPSK  SOI. 

If  the  ideal  spectral  correlation  magnitude  |S''(  /)|  used 
in  (75)-(77)  is  not  known  (to  within  a  factor  that  is  in¬ 
dependent  of/),  then  it  can  be  replaced  with  the  estimate 
(cf.  (44b))  in  (75)  to  obtain 

d"'  (T)  =  </;(T)  0  d'jT)  (78) 


where  d'„{T)  is  an  estimate  of  (74)  and  d'„(T)  is  the  esti¬ 
mate 


f 

Re  1 1  |Sk(/)|-k"'" 


d'jT) 


(  -  //j. 

|r|  <  l/2a*. 

^0,  kl  >  \/2cc,.  (79) 


F.  CLP  Method 

Since  both  SPECCORR  and  SPECCOA  methods  can 
be  interpreted  in  terms  of  indirectly  estimating  the  slope 
D  of  the  linear  phase-versus-frequency  characteristic  dis¬ 
played  by  their  integrands,  a  natural  alternative  is  to  di¬ 
rectly  estimate  this  slope  by  doing  a  linear  least  squares 
fit  to  the  phase  data  of  either  the  ratio  or  conjugate  product 
of  S“  (/)  and  5“(/).  In  both  cases,  this  phase 

('5"  <f)') 

nf)  ^  arg  =  arg  {S?.(/)S“(/)* }  (80) 

is  given  by  (cf.  (44a),  (44b)) 

'!>(/)  =  <A  -  27r(/ -I-  a/2)D.  (81) 

Thus,  we  want  to  fit  the  line  0  —  2ir{  f  -t-  u/2)D  (with 
intercept  /  =  0/2irD  —  a  12  and  slope  — 2irD)  to  the 
measured  phase  data  4>(/).  However,  without  prior 
knowledge  of  |S"(/)|,  we  do  not  know  which  band(s)  of 
frequencies  over  which  to  make  this  fit.  Nevertheless, 
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since  the  spectral  autocoherence  magnitude 


[SUf+  ol/2)  +  S“(/+  «/2)l[5?(/-  a/2)  +  5«(/-  a/2)l 


(82) 


is  close  to  unity  (cf.  Appendix)  in  bands  that  are  relatively 
clean  (S°(f  ±  a/2)  »  S'^(f  ±  a/2))  and  is  much 
smaller  than  unity  in  highly  corrupted  bands  (S^(f  ± 
a/2)  «  SU(/ ±  a/2)),  then  we  can  use  an  estimate  of 
this  function  as  a  weighting  function  W(f)  =  |C“(/)|^ 
to  perform  a  weighted  least  squares  fit  to  the  phase  data: 

min  I  \i{  f)  -  ^  +  2v(f  +  a/2)T|'IV(/)  df]^. 

(83) 

The  solution  to  this  cyclic  linear  phase  (CLP)  method  is 
given  by 

( ?(/)i7]  mf)df 

-I  J 

D  =  - - ; -  (84a) 

2w  f 

j  [  'f]~mf)df 

where 

j  4>if)W(f)df 

I  =  !>(/) - (84b) 

j  mf)df 
\fmf)df 

- •  (84c) 

]  mf)df 

In  some  cases,  the  performance  of  the  CLP  method  can 
be  improved  by  discarding  outliers  in  the  phase  data 
4>(  /)  that  are  discovered  after  the  fit  has  been  made,  and 
then  refitting  by  reusing  (84)  with  the  expurgated  phase 
data.  Also,  the  numerator  alone  in  (27)  can,  in  some  cases, 
be  a  better  weighting  function  W{  f)  since  the  corruption 
in  an  estimate  of  (80)  that  the  denominator  of  (27)  atten¬ 
uates  eventually  becomes  negligible  as  the  data  collection 
time  grows.’ 

G.  CPD  Method 

All  methods  of  TDOA  estimation  presented  up  to  this 
point  are  based  primarily  on  the  cross  correlation  or  cross¬ 
spectrum  properties  (43a)  or  (44a).  But  the  autocorrela¬ 
tion  and  autospectrum  properties  (43b)-(43c)  and  (44b)- 
(44c)  also  can  be  used  as  a  basis  for  TDOA  estimation 
since  the  phase  difference  for  either  the  two  cyclic  auto¬ 
correlations  (43b)  and  (43c)  or  the  two  cyclic  autospectral 
densities  (44b)  and  (44c)  is  simply  =  2iraD.  Thus,  as 

'The  CLP  method  without  weighting  was  first  propt>sed  by  Chen  in  |2 1 1, 
and  the  development  and  reduction  to  practice  including  the  incorporation 
of  weighting  functions,  was  carried  out  jointly  by  Gardner.  Chen,  and  M. 
F.  Kahn 


long  as  this  phase  difference  does  not  exceed  t  in  mag¬ 
nitude.  which  requires  that  the  TDOA  D,  not  exceed  1  /2a 
in  magnitude,  then  an  estimate  D  of  that  TDOA  can  be 
obtained  from  an  estimate  of  the  phase  difference  sim¬ 
ply  by  dividing  by  2ia.  In  many  practical  applications, 
for  example,  those  where  the  two  antenna  elements  are  on 
two  separate  platforms.  |D|  will  greatly  exceed  l/2a. 
Since  we  can  only  determine  D  modulo  1/a,  then  this 
phase-difference  method  will  result  in  an  ambiguity  in  the 
TDOA  estimate  except  in  applications  where  the  two  an¬ 
tenna  elements  are  on  a  single  platform  and  are  suffi¬ 
ciently  close  together  that  |D|  <  1  /2a.  However,  in  some 
situations  such  ambiguities  can  be  removed  by  other  tech¬ 
niques.  e.g.,  using  frequency  difference  of  arrival  for 
moving  platforms  or  moving  signal  sources,  or  using 
TDOA  at  different  times  for  moving  platforms. 

Paralleling  the  derivation  (65)-(66)  of  the  SPECCOA 
method,  we  can  exploit  the  phase-difference  property 

RUu)  =  C-/?“(M)c"'-"“^  (85) 

where  C‘  =  |/4|',  over  a  range  of  lag  values  u  by  seeking 
the  value  of  the  estimates  t  =  D  and  C  that  minimize  the 
sum  of  squares  of  error  magnitudes  between  the  measure¬ 
ments  of  the  left  and  right  sides  of  (85)  with  D  replaced 
by  t: 

D  =  arg  min  {^(t)}  (86) 

i'.r 

where  e„(T)  is  an  estimate  of 

ejr)  =  J  |/??(M)  -  C‘R'/{u)e''‘^"^\-  du.  (87) 

The  result  is  given  explicitly  by 

D  =  ^  angle  |  j  /??(«)/?“(«)*  rfuj  (88a) 

=  ^  angle  J  5“  ( /)  5“  ( /)*  j  (88b) 

and  is  referred  to  as  the  cyclic  phase-difference  (CPD) 
method. 

Unlike  the  CCCC,  SPECCORR,  SPECCOA.  SPEC- 
CON,  and  CLP  methods  preceding  this  CPD  method,  all 
of  which  provide  genuine  time-differer  ;e  estimates,  the 
CPD  method  provides  phase-difference  estimates  only. 
That  is,  it  regenerates  sine  waves  at  each  of  the  sensors 
(by  forming  the  lag  products  x(t  +  u/2)x(t  —  u/2)  and 
,v(/  -t-  u/2)y(t  -  «/2)),  then  measures  their  magnitudes 
and  phases  (by  averaging  the  product  of  each  of  the  lag 

"The  CPD  melhixi  was  invented  by  Gardner  and  first  disclosed  in  a  1989 
report  for  a  classified  research  project.  It  was  reduced  to  practice  jointly 
by  Gardner  and  C.  M.  Spooner  and  is  being  publicly  disclosed  for  the  first 
lime  in  this  paper. 
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products  with  the  complex  sine  wave  exp  (-ilirat)  to  ob¬ 
tain  /?v(m)  and  /?“(w)).  then  averages  the  conjugate  prod¬ 
ucts  Ry(u)Rj{u)*  of  these  complex  quantities,  each  prod¬ 
uct  of  which  ideally  has  the  same  phase  IvaD,  and, 
finally,  extracts  the  resultant  phase  of  this  average.  In  fact, 
it  can  be  shown  that  the  CPD  method  can  be  interpreted 
as  a  phase  interferometiy  method  since  the  same  estimate 
(88)  results  from  nulling  out  the  regenerated  sine  wave  in 
the  difference  of  lag  products  by  minimizing  the  sum  of 
squared  magnitudes  of  the  Fourier  coefficients  of  the  lag- 
product  difference  with  respect  to  an  adjustable  time-dif¬ 
ference  (or.  effectively,  phase-difference)  parameter  t: 

D  =  arg  min  {/„(t)}  (89) 

T 

where 

/.(’■)=  I  I  <  [C'.xU  ~  T  +  u/2)xu  -  T  -  ujl) 

-  yu  -t-  M/2).v(f  -  u/2)le  >  |-  du  (90) 

in  which  <  •  >  is  a  finite-time  average.  In  comparison,  the 
SPECCON  method  nulls  the  sum  of  squared  magnitudes 
of  the  Fourier  coefficients  (7 1 )  of  the  lag  products  of  the 
difference  signals  x{t)  -  y(/  +  t)  with  respect  to  an  ad¬ 
justable  time-difference  parameter  r.  Thus,  the  SPEC- 
CON  method  forms  difference  signals,  one  member  of 
which  contains  an  adjustable  delay  parameter,  before  sine- 
wave  regeneration,  whereas  the  CPD  method  forms  dif¬ 
ference  signals,  one  member  of  which  contains  an  adjust¬ 
able  delay,  after  sine-wave  regeneration. 

H.  Generalization  to  More  Than  Two  Receivers" 

When  more  than  two  receivers  corresponding  to  mul¬ 
tiple  reception  platforms  and/or  multiple  antenna  ele¬ 
ments  on  a  single  platform  are  available,  each  distinct  pair 
of  receivers  can  be  used  to  obtain  a  distinct  estimate  of 
TDOA.  If  there  are  N  receivers,  then  there  are  M  =  (N 
-  l)N/2  distinct  receiver  pairs.  If  the  distance  between 
the  antenna  elements  for  the  mth  pair  is  L„  and  the  broad¬ 
side  angle  of  the  /nth  pair  with  respect  to  an  arbitrarily 
given  direction  is  d„,  then  the  TDOA  D„  resulting  from 
a  far-field  signal  with  a  linear  wavefront  arriving  from 
angle  $  relative  to  the  same  given  direction  is 

D„  =  —  sin  (6  -  e„).  (91) 

c 

Consequently,  the  Af  TDOA  estimates  D„  obtained  by  any 
method(s)  can  be  combined  to  obtain  a  single  AOA  esti¬ 
mate  d  as  follows: 

where  <  •  >  denotes  average  over  m  =  1,  2,  •  •  •  ,  M.  Of 

"am  TDOA-based  methods  for  .AOA  estimation  using  multiple  pairs  of 
sensors  !ha!  are  presenied  in  fhi.s  .sub.section  were  invenled  by  Gardner  and 
reduced  to  practice  jointly  by  Gardner  and  C.  M.  Spooner  and  are  being 
publicly  disclosed  for  the  first  time  in  this  paper. 
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course,  if  some  of  the  M  individual  AOA  estimates  are 
expected  to  be  more  accurate  (e  g.,  because  the  corre¬ 
sponding  distance  between  antenna  elements  is  larger), 
then  these  can  be  weighted  more  heavily  in  a  weighted 
average  over  m. 

In  order  to  remove  the  ambiguities  (due  to  ±)  in  (92), 
a  clustering  technique  can  be  applied  to  the  set  of  2M  es¬ 
timates  produced  by  (92)  to  select  the  appropriate  M  es¬ 
timates  (i.e.,  to  remove  the  ambiguities)  prior  to  averag¬ 
ing.  It  should  be  noted  that  when  some  sensor  pairs  have 
especially  large  separations,  all  other  pairs  should  prob¬ 
ably  be  omitted  since  their  TDOA  estimates  will  (all  else 
being  equal)  be  considerably  less  accurate.  When  the  far- 
field  assumption  of  linear  wavefronts  is  violated,  as  it 
often  is  for  widely  separated  platforms,  the  following  al¬ 
ternative  method  can  be  used.  This  alternative  also  avoids 
the  need  for  clustenng  to  remove  ambiguity. 

Although  the  objective  in  this  paper  is  to  introduce  new 
methods  for  TDOA  estimation,  which  is  usually  viewed 
as  an  alternative  to  the  high-resolution  direction-finding 
approach  based  on  antenna  arrays  with  closely  spaced  ele¬ 
ments  (e.g.,  fraction-of-a-wavelength  spacing),  it  .should 
be  mentioned  that  there  is  some  potential  for  deriving  new 
array-based  methods  from  TDOA  methods.  This  is  espe¬ 
cially  evident  for  the  CPD  method  since  it  is  actually  a 
phase-difference  method.  In  fact,  the  least  squares  ap¬ 
proach  that  led  us  to  the  CPD  algorithm  can  be  general¬ 
ized  from  A/  =  2  to  A  >  2.  and  the  resultant  algorithm'" 
shows  some  promise  for  high-resolution  signal-selective 
source  location.  But.  as  explained  above,  the  CPD  method 
is,  like  other  array-based  methods,  restricted  to  closely 
spaced  sensors  if  ambiguities  are  to  be  avoided.  In  con¬ 
trast  to  this,  some  of  the  TDOA  methods  described  in  this 
paper  can  be  generalized  from  sensor  pairs  to  .sensor  ar¬ 
rays  (multiple  sensor  pairs  or  platform  pairs)  and  there  is 
no  restriction  on  sensor  separation  for  these  methods. 
Moreover,  unlike  typical  array-based  methods,  which  are 
designed  for  narrow-band  signals,  the  TDOA-based 
methods  apply  directly  to  wide-band  signals  as  well  as 
narrow-band  signals.  There  are  two  approaches  to  accom¬ 
plishing  this  generalization.  The  first  approach  is  to  pro¬ 
duce  separate  TDOA  estimates  and  then  combine  them  to 
obtain  a  single  AOA  estimate  as  in  (92).  The  second  ap¬ 
proach  combines  cyclic  correlation  or  cyclic  spectrum 
measurements  for  multiple  sensor  pairs  and  uses  these 
jointly  and  directly  to  estimate  a  single  AOA.  Since  all 
measurements  are  used  jointly  in  a  single  optimization  in 
the  second  approach,  whereas  subsets  of  measurements 
are  used  individually  in  each  of  multiple  optimizations  in 
the  first  approach,  the  second  approach  is  expected  to  yield 
superior  methods.  Nevertheless,  preliminary  simulations 
suggest  that  even  these  superior  TDOA-based  methods 
will  be  inferior  to  high-resolution  direction-finding  meth¬ 
ods  when  they  are  applicable,  that  is,  when  the  sensors 
are  separated  by  only  fractions  of  a  wavelength. 

'■ThLs  cyclic-pha.se  algorithm  wa.s  hrsi  de.scribed  by  .Xu  and  Kailalh  |22|. 
and  is  put  in  perspective  relative  to  other  cyclo.stationarily-cxploiting  array- 
based  methods  in  |7), 
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As  one  example  of  the  multiple  sensor-pair  methods, 
the  least  squares  SPECCORR  method  generalizes  from 
(61)  to  the  M-SPECCORR  method'’ 

d  =  arg  max  b'^'  sin  [6  “  j  (93) 

where  b^'  (t)  is  an  estimate  of  b'^{T)  for  the  mth  pair  of 
sensors  which  produce  the  received  signals  jc„(f)  and  >„,(/) 
in  place  of  T(f)  and  v(r)  in  (62). 

Similarly,  thf  least  squares  SPECCOA  method  gener¬ 
alizes  from  (66  :o  the  M-SPECCOA  method  (assuming 
real  data) 

0  =  arg  max  c”'  sin  10  “  j  (94) 

where  is  an  estimate  of  c^(t)  for  the  wth  pair  of 

sensors. 

When  the  separation  between  sensors  is  large  (e.g., 
when  the  sensors  are  on  separate  platforms),  the  assump¬ 
tion  of  plane-wave  (far-field)  reception  can  break  down, 
and  the  unique  relationship  (91)  between  ±AOA  and 
TDOA  can  fail  to  hold.  In  this  case  the  hyperbolic  rela¬ 
tionship  between  all  source  locations  (AOA  and  range) 
that  give  rise  to  a  specific  TDOA  must  be  used  in  place 
of  (91).  If  the  two  sensors  in  the  mth  pair  have  ranges  of 
Ti  and  and  angles  of  0|  and  01  relative  to  an  arbitrarily 
given  reference  point,  and  if  the  source  has  range  r  and 
angle  0,  relative  to  this  same  reference  point  then  we  have 

D„  =  ±  -  iKr,  cos  0|  -  r  cos  0)’ 
c 

+  (r,  sin  0|  -  r  sin  0)^1 

-  [(r2  cos  01  -  r  cos  0)^ 

+  (r2  sin  01  -  r  sin  0)^]'^^j  (95) 

where  c  is  the  speed  of  propagation  and  D„  is  the  TDOA, 
The  loci  of  all  range/angle  pairs  (r,  0)  that  satisfy  (95)  for 
each  TDOA  D„  is  a  hyperbola.  If  the  ±  sign  in  (95)  is 
known  by  virtue  of  the  orientation  of  the  sensor  pair,  then 
we  can  simply  substitute  (95)  with  the  appropriate  sign  in 
place  of  D„  =  (L„Jc)  sin  (0  —  0„)  in  the  M-SPECCORR 
or  M-SPECCOA  algorithms  and  replace  the  univariate 
search  over  0  with  a  bivariate  search  over  r  and  0,  and 
thereby  estimate  both  range  and  AOA.  It  is  pointed  out 
that  this  joint  optimization  approach  to  emitter  location 
should  be  contrasted  with  the  conventional  approach  in 
which  individual  TDOA’s  are  estimated  and  the  corre¬ 
sponding  hyperbolic  loci  of  range/angle  pairs  are  inter¬ 
sected  to  remove  ambiguity  and  determine  a  unique  em¬ 
itter  location. 

In  closing  this  section,  it  is  pointed  out  that  the  perfor¬ 
mance  of  all  TDOA  methods  presented  here  can  be  de¬ 
graded  in  the  unusual  situation  where  more  than  one  spec¬ 
trally  overlapping  signal  with  the  same  cycle  frequency  a 

'  ‘As  an  alternative,  the  .sum  in  (9,t)  can  be  replaced  with  a  product. 


is  received.  (This  occurs,  for  example  in  a  multiuser 
code-division-multiplexed  communication  system.)  How¬ 
ever.  the  CCC  method  based  on  (49),  CCCC  method 
based  on  (51)  (and  that  based  on  (54)),  and  the  SPECCOA 
method  based  on  (67)  can  operate  properly  when  there  are 
multiple  independent  identically  distributed  (except  for  a 
complex  scale  factor  and  time  of  arrival)  spectrally  over¬ 
lapping  SOI  with  the  same  a  provided  that  their  TDOAs 
are  resolvable  (cf.  Section  I).  Moreover,  the  CCC  method 
can  operate  properly  even  if  the  multiple  signals  are  not 
identically  distributed.  Also,  the  generalization  of  the 
CPD  method  to  N  >  2  sensors  and  the  other  cyclosta- 
tionarily -exploiting  array-based  methods  [7]  can  operate 
properly  provided  that  the  number  of  signals  with  cycle 
frequency  a  is  less  than  N. 

V.  Conclusion 

We  have  introduced  a  new  class  of  methods  for  signal- 
selective  TDOA  estimation.  By  exploiting  the  cyclosta- 
tionarity  property  of  the  signal  of  interest,  as  reflected  in 
the  spectral  correlation  functions  for  the  received  data, 
the  effects  of  additive  noise  and  interfering  signals  are 
ideally  (for  unlimited  data  collection  times)  removed  by 
these  methods.  As  shown  in  part  II  (1)  of  this  two-part 
paper,  these  methods  can  be  highly  tolerant  to  noise  and 
interference  for  finite  data-collection  times.  The  signal  se¬ 
lectivity  of  these  methods  makes  it  possible  to  apply  them 
to  closely  space  as  well  as  widely  separated  sensors  and 
to  generalize  them  to  multiple  sensor  pairs,  including 
multiple  platforms  and  possibly  arrays  of  closely  spaced 
sensors.  Since  their  principles  of  operation  are  indepen¬ 
dent  of  signal  bandwidth  (assuming  that  the  propagation 
channel  and  sensor  characteristics  are  either  matched  or 
are  independent  of  frequency  throughout  the  signal  band), 
these  methods  apply  to  both  narrow-band  and  wide-band 
signals.  Given  this  general  applicability,  it  can  be  seen 
that  the  new  class  of  methods  at  least  partially  bridges  the 
gap  between  the  previously  dichotomous  array-based  di¬ 
rection  finding  problem  and  the  time-difference  location 
problem. 

Appendix 

The  BPSK  Signal''* 

A  BPSK  signal  can  be  modeled  by 

s{t)  =  a(t)  cos  (2^/7  +  <i>o)  (Al) 

where 

00 

a{t)  =  E  a„q(t  -  to  -  nT^)  (A2) 

n  =  -  00 

in  which  q(t)  is  a  finite-energy  keying  envelope  with  peak 
height  of  unity  and  width  of  7'<.=  l/a*.  and  {«„ }  is  a 
binary  (±1)  sequence.  If  {cr„ }  is  modeled  as  a  sequence 
of  independent  variables  with  equiprobable  values,  then 

'■’This  Appendix  is  an  abbrcvialion  of  the  more  complete  treatment  given 
in  (10). 
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it  can  be  shown  that  the  spectral  correlation  function  is 
given  by 

5:(/)  =  T^r  le  f  /  +  /o  +  e*  f  /  +  /o  -  ^ 


4^  1 


a 


a 


for  a  =  n  /  and 


,  -j2Taio 


s?(/) 


1 

47-* 


.  -i(2i|q  2/;  Id)  2*11 


+  47  (/-/,  +  ^  (2*  (/  +  /,-  “ 


•  e 


-f(2irla  -  2/i  |m  -  2<&o> 


(A3) 


fora  =  ±2/  +  n/T’ji  for  all  integers  n.  and  S'^(f)  =  0 
for  all  other  values  of  a.  In  (A3),  Q(f)  is  the  keying 
envelope  transform 


-  iZtfl 


dt. 


SOD 

q(t)e 

—  OD 

The  inverse  Fourier  transform  of  (A3)  is  given  by 

JQO 

df 

—  00 


(A4) 


=  r"(r)  cos  (2t/.t)c 

lit 


lope,  but  those  in  Figs.  2-4  were  obtained  using  the  more 
practical  full-duty-cycle  half-cosine  keying  envelope: 


q(t)  = 


I'l  ^  TJ2 

otherwise. 


(A8) 
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where 

/■"(r) 


for  a  =  n/T*  for  all  integers  n. 

Graphs  of  the  magnitudes  of  the  functions  S“(/)  and 
/?“(t)  are  shown  in  Figs.  5  and  6  for  the  case  of  a  full- 
duty-cycle  rectangular  keying  envelope 


qU)  ^ 


1^1  ^  TJ2 
otherwise 


(A7) 


with  keying  rate  a^  =  1  /  =  1/10  and  carrier  frequency 

/r  =  0.33.  In  these  graphs  the  magnitudes  of  the  complex¬ 
valued  functions  are  plotted  as  the  heights  of  surfaces 
above  the  planes  with  coordinates /and  a.  Since  the  func¬ 
tions  can  be  nonzero  for  only  discrete  values  of  the  cycle 
frequency  parameter  a,  the  surfaces  consist  of  infinitesi¬ 
mally  thin  slices  parallel  to  the /or  t  axis. 

The  graphs  of  the  TDOA  functions  for  BPSK  shown  in 
Fig.  1  were  obtained  using  the  rectangular  keying  enve- 
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